##################################################################

# Gautam Nair
# gautam.nair@yale.edu
# Misperceptions of Relative Affluence and Support for International Transfers
# Make Table 10: Observational Correlation between Preferences and Misperceptions

##################################################################

# setwd("")

##################################################################
# loading packages
##################################################################

library(Hmisc)
library(foreign)
library(sandwich)
library(lmtest)
library(numDeriv)
library(stargazer)
library(ggplot2)
library(plyr)
library(gridExtra)
library(ri)
library(dplyr)
library(plyr)
library(scales)

##################################################################

rm(list=ls())

##################################################################
#Observational Correlation between Preferences and Misperceptions
##################################################################

data.working <- readRDS("d_r_cleaned_analysis_dataset.Rda")

data.working <- data.working[data.working$group2==1| data.working$group3==1,]

#data.working <- data.working[data.working$group2==1,]

data.working$rank.estimate.diff.10 <- data.working$rank.estimate.diff/10
data.working$median.estimate.top.100k.10000 <- data.working$median.estimate.top.100k.10000 - ((2100/10000))
data.working$median.estimate.top.100k.10000

vars1 <- c(
"rank.estimate.diff.10")

vars2 <- c(
"rank.estimate.diff.10",
"charity.pretreat.yes",
"income.new",
"age",
"education.categorical",
"ideo.new",
"identity.cosmopolitan",
"partyid",
"religion.attendance",
"work.employed",
"gender.female",
"group3"
)

vars3 <- c(
"median.estimate.top.100k.10000")

vars4 <- c(
"median.estimate.top.100k.10000",
"charity.pretreat.yes",
"income.new",
"age",
"education.categorical",
"ideo.new",
"identity.cosmopolitan",
"partyid",
"religion.attendance",
"work.employed",
"gender.female",
"group3"
)

varsomit <- c(
"income.new",
"charity.pretreat.yes",
"age",
"education.categorical",
"ideo.new",
"identity.cosmopolitan",
"partyid",
"religion.attendance",
"work.employed",
"gender.female",
"group3"
)

list.ind.group <- list(vars1, vars2, vars3, vars4)

dep.var <- c("q13ab.3", "q8.toolittle", "q11.yes")
reg.models <-vector("list", length(dep.var)) 
se.models <-vector("list", length(dep.var))
matrix <- vector("list", length(dep.var))

x=0 

for(i in 1:length(dep.var)){
	for(y in 1:length(list.ind.group)){
		my.formula <-paste(dep.var[i],'~', paste(list.ind.group[[y]], collapse= ' + '))
		temp.model <- lm(my.formula, data=data.working)
		temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
		x = x+1
		se.models[[x]] <- temp.se
		reg.models[[x]] <- temp.model
	}
}

varlabels <-       c("Bias: Rank", "Bias: Median")
spectitle <- c("Correlation between Misperceptions of Relative Income and Preferences for International Redistribution")
outputfile <- c("tf_t_10_misperceptions_preferences_correlation.txt")
columnlabels <-c("Donation to Int. Charity (%)", "Increase Foreign Aid (1/0)", "Reduce Ag. Trade Protections (1/0)")
sepcolumns <- c(4,4,4)
cutoffs.star <- c(0.05, 0.01)
char.star = c("**", "***")
lines.add <- c("Controls", "N", "Y", "N", "Y","N", "Y","N", "Y","N", "Y","N", "Y")
specnotes <- c("Robust standard errors in parentheses",
          "Dependent variable in (1) takes values of:",
           "1 (not a problem at all)",
           "2 (a small problem)",
           "3 (a problem)",
           "4 (a serious problem)",
           "5 (a very serious problem)",
          "Dependent variable in (2) takes values of:",
          "1 (serious or very serious problem), 0 (otherwise)",     
           "Control group means = (1) 3.151 (2) 0.394")

stargazer(reg.models[[1]],  reg.models[[2]], reg.models[[3]], reg.models[[4]], reg.models[[5]], reg.models[[6]], reg.models[[7]], reg.models[[8]],
reg.models[[9]], reg.models[[10]], reg.models[[11]], reg.models[[12]],
          se=list(se.models[[1]],  se.models[[2]], se.models[[3]], se.models[[4]], se.models[[5]], se.models[[6]], se.models[[7]], se.models[[8]], se.models[[9]], se.models[[10]], se.models[[11]], se.models[[12]]),
         title= spectitle,
	 type="text",
          out=outputfile,
          no.space=TRUE, 
          model.numbers= TRUE,
          align=FALSE,
        covariate.labels=varlabels,
          header= FALSE,
          font.size="small",
          model.names=TRUE,
          digits=2,
          omit.stat = c("rsq", "f","adj.rsq"),
         column.labels  = columnlabels,
         column.separate =sepcolumns,
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
         df=FALSE,
         omit.table.layout="n",
                  star.cutoffs = cutoffs.star,
     star.char = char.star,
	add.lines=list(lines.add),
	column.sep.width = "0.5pt",
         omit= varsomit,
         notes.align=c("l"),
notes.label = "",
 notes        = specnotes,
  notes.append = FALSE

)   
 
